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ABSTRACT 


Motion  control  of  underway  replenishment  operations  is 
achieved  through  the  use  of  sliding  mode  control  with  a  Linear 
Quadratic  Gaussian  compensator  design.  External  disturbances 
include  first-order  wave  force  and  moment  as  well  as  slowly 
varying  interaction  forces  and  moments  between  the  two  ships. 
Feedback  control  is  used  to  provide  adequate  stability  of 
motions  while  feedforward  control  with  disturbance  estimation 
and  compensation  achieves  the  desired  steady  state  accuracy. 
The  results  demonstrate  that  satisfactory  path  keeping  during 
operations  can  be  maintained  for  various  ship  proximity 
distances  and  environmental  conditions. 
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I .  INTRODUCTION 


A.  AIM  OF  THIS  STUDY 


Accurate  course  keeping  is  one  of  the  most  important 
problems  for  ship  maneuvering  and  especially  in  an  underway 
replenishment  in  an  open  sea.  There  is  an  apparf it  danger  in 
such  an  operation  due  to  the  complex  effects  of  the  seaway, 
and  the  conventionally  slow  control  loop  indicated  in  Figure 
1.1. 


Figure  1.1  Close  Loop  Control  System  for  Ship  Maneuvering 


In  particular,  during  underway  replenishment  the  ship  is 
subject  to  a  variety  of  forces  and  environmental  disturbances. 
The  most  important  of  these  are  the  interaction  forces  and 
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moments  that  are  generated  as  a  result  of  the  two  ships  in 
proximity  with  one  another [Ref . 7] [Ref. 8] .  These  forces  depend 
on  the  relative  positions  of  the  two  ships  and  both  their 
magnitude  and  sign  is  not  easily  predictable  by  analytical 
techniques  In  a  random  seaway  there  exist  also  first-order 
irregular  wave  forces  and  second-order  wave  drift  forces.  The 
latter  appear  in  the  form  of  long  term  sT  owly  varying 
excitations  and  do  not  have  a  significant  effect  on  lateral 
separation  control  [Re^.' .  8]  .  In  this  thesis  we  focus  our 
attention  on  the  effects  of  the  interaction  forces  and  first- 
order  wave  forces  on  the  system  response.  We  employ  a  sliding 
mode  feedback  control  law  due  to  its  robustness 
characteristics  with  regard  to  unmodeled  dynamics  and 
disturbances.  The  sliding  surface  is  based  on  the  Linear 
Quadratic  Regulator.  Partial  staee  measurement  is  assumed  and 
a  Kalman  filter  is  designed  to  provide  an  accurate  estimate  of 
the  unmeasurable  states  and  to  minimize  the  effects  of 
measurement  noises.  Disturbance  estimation  and  compensation  is 
used  in  order  to  ensure  the  desired  steady  state  accuracy.  A 
schematic  block  diagram  of  the  control  design  is  shown  in 
Figure  1.2. 
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1.2  Automatic.  Contrc."'  of  Ship  Replenishment  at  Sea 


3 


B.  THESIS  OUTLINE 


In  Chapter  II,  the  LQR  method  of  sliding  mode  controller 
design  is  developed  and  simulated  at  calm  sea.  Chapter  III  is 
conserned  with  the  Kalman  filter  design.  Simulations  of  the 
UNREP  with  open  sea  effects  are  presented  in  Chapter  III. 
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II.  MATHEMATICAL  MODEL  IK  THE  HORIZONTAL  PLANE 


A.  NONLINEAR  STEERING  EQUATIONS  IN  THE  OPEN  SEA 

The  UNREP  operation  at  sea  is  schematically  depicted  in 
Figure  2.1  where  Y  and  N  are  the  interaction  forces  and 
moments.  The  leading  ship  keeps  its  course  unchanged  and  the 
tracking  ship  keeps  the  separation  distance  constant,  normally 
less  then  100  feet  (depending  on  the  type  of  ship) . 

Figure  2.6  and  Figure  2.7  show  the  interaction  forces  and 
moments  on  the  tracking  ship  due  to  the  interaction  effect  of 
the  leading  ship.  The  interaction  force  Y(A,B)  acts  through 
the  ship's  center  of  gravity  and  the  moment  N(A,B)  about  the 
center  of  gravity  in  the  horizontal  plane.  The  two  curves  used 
in  this  study  are  for  a  typical  Mariner  ship  at  a  nominal 
speed  of  15  knots  and  for  different  relative  positions.  In 
applying  these  curves  to  the  leading  ship,  the  interaction 
force  and  moment  need  to  be  changed  to  -Y(-A,B)  and  -N(-A,B), 
respectively.  For  each  figure,  the  curves  for  B=50  ft  and 
B=100  ft  were  determined  from  experimental  model  testing 
results [Ref . 6] [Ref. 8],  and  the  curves  in  between  were 
determined  by  interpolation.  The  data  were  also  extrapolated 
up  to  B=150  ft.  A  and  B  are  continually  read  by  the  simulation 
program  as  input  data  to  a  two  dimensional  interpolation 
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routine  which  calculates  the  base  line  interaction  forces  and 
moments  acting  on  the  ship. 

The  coordinate  system  used  in  this  study  follows  the 
standard  system  from  Principles  of  Naval  Architecture  [Ref .  11]  . 
In  practice,  the  leading  ship  keeps  constant  heading,  while 
the  tracking  ship  maintains  a  constant  lateral  separation 
distance.  Different  configurations  and  the  directions  of 
variations  of  the  interaction  force  and  moment  as  the  tracking 
ship  moves  from  -524  ft  to  524  ft  are  shown  in  Figure  2.2. 

The  primary  consern  for  UNREP  simulation  is  related  to  the 
space  coordinates  (longitudinal  and  lateral  separation 
distance,  yaw  angle,  and  yaw  rate)  .  So  the  equations  of  motion 
are  first  expressed  in  ship  velocity  coordinates  (u,  v)  and 
then  transformed  to  space  coordinates (x,  y)  (Figure  2.3).  The 
transformation  equations  are 

x=ucosi|f-vsini|f  (2.1) 

y=usintjr  + vcosil;  (2.2) 

The  ship's  coordinates  are  assumed  to  be  at  the  ship's 
center  of  gravity  (Xg=0,  yc=0  )  . 
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Passes  the  Leading  Ship 


Note  : 

A:  longitudinal  separation  distance. 

E:  lateral  separation  distance. 

Arrows  indicate  positive  directions  according  to 
established  sign  conventionIRef.il]. 
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Figure  2,2  Relative  Positions  of  The  Two  Ships  for  A=-524 
ft,  0,  and  +524  ft  and  Corresponding  Interaction 
Forces  and  Moments . 

*  Note:  LS :  Leading  Ship. 

TS :  Tracking  Ship. 
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POSITION  OF  CENTER  Of  GRAVITY 
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B.  MANEUVERING  MATHEMATICAL  MODEL 

The  nonlinear  maneuvering  equations  of  motion  and  the 
hydrodynamic  interaction  forces  and  moments  that  act  on  the 
ship's  hull  are  modeled  for  15  knots  nondimens ional  forward 
speed.  In  an  UNREP  control  simulation,  the  ability  to  manually 
or  automatically  control  the  rudder  and  propeller  speed  must 
be  present.  Only  rudder  input  is  considered  here,  and  the  ship 
forward  speed  is  assumed  constant  by  the  propulsion  control. 
The  basic  mariner  study  ship  characteristics  are  presented  at 
Table  2.1. 

Table  2.1  Characteristics  of  Mariner-type  Study  Ship 


Length  between  perpendicular 

527.8  ft (160.87  m) 

Beam 

76.0  ft(  23.26  m) 

Draft 

29.75  ft (  9.06  m) 

Displacement 

16800.0  tons 
(15120  t) 

Block  coef  f  icient  (Cf,) 

0.6 

*  The  ship's  coordinate  system  is  assumed  to  be  at  the 


ship's  center  of  gravity(xG,  yc=  0)  . 
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The  goal  of  this  study  is  to  present  an  alternative  way  of 
steering  control  based  on  sliding  modes.  For  this  reason  the 
simulation  does  not  concentrate  on  the  derivation  of  ship 
dynamic  and  hydrodynamic  effects,  but  on  the  results  of 
control  and  observer  feedback.  Only  the  tracking  ship  is  the 
major  consideration,  and  we  assume  a  reliable  lateral  distance 
measurement.  The  ship  hydrodynamic  maneuvering  model  for  the 
study  ship  consists  of  the  nonlinear  equations  of  motion  in 
the  horizontal  plane  (sway  and  yaw) .  The  hydrodynamic 
coefficients  are  presented  in  Table  2.3,  and  the 
nondimens ional  equations  of  motion  are: 

SWAY  EQUATION: 

v+  imx^-Y^)  r=F2  (u,v,  r,6)  (2.3) 


YAW  EQUATION: 


(mx^-N^,)  +  r=F^  (u,  V,  r,  6) 


(2.4) 


where 


F^iu,  V.  r ,  b)  =Y^v+ (Y^-wu^)  r+Y^6-*-Y^p  +  Y{A,  B) +  ^2  5) 


and 
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Fj  (u,  V,  r,  6)  =N^v+  (Nj.-mxQU^)  r+Ni,6+N^n+N{A,  B)  +n^+n2  (2.6) 


Y(A,B)  nordimensional  hydrodyricimic  interaction  force 
caused  by  leading  ship. 

N(A,B)  nondimens ional  hydrodynamic  interaction  moment 
caused  by  leading  ship. 
f,,n,  first-order  component  of  random  sea. 
f2,n2  second- order  component  of  random  sea. 

Some  assumptions  were  made  for  simplification  without  any 
significant  physical  effects.  First,  the  ship's  coordinates 
are  assumed  to  be  at  the  ship's  center  of  gravity  {  x^,  « 

0)  ;  second,  Y„n  and  N„n  are  small  terms  also  neglected  from  the 
simulation;  and  third,  the  second-order  force  and  sway  moment 
£2  and  nj  are  not  considered  in  this  thesis.  This  is  because 
the  effect  of  these  slowly  varying  forces  and  moments  is 
similar  to  the  interaction  forces  and  moments  which  are 
analyzed  in  detail. 

1.  Rudder  Dynamics  Model 

There  are  two  limiters  for  modeling  rudder  dynamics, 
the  first  limiter  models  the  rudder  saturation  limits  at  ±0.4 
radians.  The  second  limiter  models  the  proportional  band  of  a 
variable  -  displacement  pump  by  limiting  its  maximum  percent 
stroke,  the  limits  for  this  nonlinear  element  have  been  found 
to  be  0.297  radians  for  this  study  case. [Ref. 1] [Ref. 6] 
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The  lag  time  and  "dead  band"  in  the  rudder  dynamics 
when  a  helm  angle  is  commanded  are  important  aspects  of  the 
manuvering  control  problem.  The  rate  of  change  of  the  rudder 
angle  was  assumed  to  be  inversely  proportional  to  the  error 
signal,  subject  to  a  maximum  rudder  rate  of  .0698  rad/sec,  and 
a  maximum  error  input  of  0.297  radians. 

The  system  gain  can  be  defined  as: 

k,  =  Maximum  rudder  rate/Maximum  error  input 

=  0.0689/0.298  =  0.234  /sec 
We  nondimensionalize  this  system  gain  as: 
k/=k,*L/ui=  0.234*527.8  ft/ (15*1 . 689  )  ft/sec  =4.89 
Where 

L  =  ship  length=  527.8  ft 

u,=  ship  forward  speed=  15  knots=  15*1.689  ft/sec 
This  again  was  taken  approximately  equal  to  5  in  this 
study;  and  t,  =  1/k/  is  defined  as  the  nondimensional  rudder 
system  time  constant. 

The  rudder  "dead  band"  was  set  to  ±  0.5  deg.  An  analog 
diagram  of  the  rudder  dynamics  is  shown  in  Figure  2.4  [Ref. 8], 
which  includes  the  response  of  the  rudder  system  to  step 
commands.  The  equivalent  rudder  dynamics  are  described  by  a 
first  order  lag,  see  Figure  2.5,  where 
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(2.7) 


6  =  + 


7,  is  the  rudder  dynamics  time  constant,  and 
6,  is  the  commanded  rudder  angle. 

2 .  State  Space  Representation 

The  mathematical  model  consists  of  the  previous  sway 
and  yaw  equations  of  motion 


v+  {wx^-Yf.)  i  =  Y^v*  (  Y^-mu)  x-^Yf^b+Y{A,  B)  (2.8) 

imx^-N^.)  v’+  t=N^v+  (Nj.~mxQU)  r+N(^6+N{A,  B)  +n^  (2.9) 

and  the  rudder  dynamics 

6  =  -  — 6  +  — 6_  (2.10) 

-tr 

The  lateral  separation  rate  can  be  calculated  from  the 
components  of  lateral  direction  of  ship  forward  speed  (u) 
lateral  speed  (v)  ,  and  heading  angle  (\}/) 

y=usinijf  +  vcosi|;  (2.11) 
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With  the  small  angle  assumption  we  can  linearize  this 


equation  as ; 


y=i}r+v 


(2.12) 


where  the  dimensionless  forward  speed  is  one. 

The  last  kinematic  relation  is  between  the  ship 
heading  \l/  and  heading  rate  \j/=r.  The  first-order  wave  force  f, 
and  moment  n,  in  (2.8)  and  (2.9)  vary  much  rapidly  compared  to 
the  dynamics  of  the  ship  and  for  control/observer  design,  they 
can  be  effectively  modeled  as  white  noise  components.  This  is, 
however,  not  true  for  the  interaction  forces  and  moments 
Y(A,B)  and  N(A,B).  These  vary  at  approximately  the  same  rate 
as  the  ship  dynamics  since  tb.j  uepend  on  the  relative 
separation  between  the  '.wo  ships.  They  are  best  modeled  as 
colored  noise  components?  the  response  of  a  first-order 
shaping  filter  driven  by  white  noise 

T^=-N^W^  (2.13) 

TyY=-Y*Wy  (2.14) 

where  the  dimensionless  time  constants  T^  and  Ty  are  equal  to 
l(i.e.  the  time  that  it  takes  to  travel  one  ship  length). 
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Collecting  the  above  equations  we  get  a  seven  order  system  in 
state  space  form.  The  control  input  is  the  commanded  rudder 
angle  6^  and  the  augmented  state  vector  contains  the  actual 
rudder  angle  6,  heading  angle  rp,  sway  velocity  v,  yaw  rate  r, 
lateral  separation  distance  y,  and  the  two  components  of  the 
interaction  force  Y  and  moment  N.  In  a  compact  vector 
notation,  the  complete  state  space  system  is  written  as 
follows . 


x= Ax+ B  Wj^ + Bj  W2 


(2.15) 


where 


5' 


y 

N 

y 


^2.16) 


w,=l 


N 


(2.17) 
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(2.18) 


A=[a,,], 


i=l-7; 

J=l-7. 


(2.19) 


^1  = 


0 

0 

0 

0 

0 

0 


(2.20) 


Bo 


0  0 
0  0 
0  0 
0  0 
0  0 


(2.21) 


B3=i 


Ai 

Az 

Ai 

A  2 

Ai 

A  2 

Ai 

A  2 

Az 

Ai 

A  2 

Ai 

Az 

(2.22) 
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1 


(2.23) 


a 


11 


X 


r 


^24  ^ 


-N,  ( 


_ _ Y^iNj-I^)  -N^iYt-mXs) 

( Y^-m)  iNi-I^)  -  ( Y^-mX^)  (N^-mXs) 


(Y^-mu)  iN^-I^)  -  {Yf-mX^)  {N^-mX^u) 
(  Y^-m)  iN^-I^)  -  (  Y^-mX^)  {N^-mXa) 


^ _ yf-^G _ 

(Y^-m)  iN^-I^)  -(Yf -mX^)  ( N^-mX^) 


-  _ _ _ 

( Y^-m)  (N^-I^)-{  Yf-mX^)  {N^~mXj) 


(2.24) 


(2.25) 


(2.26) 


(2.27) 


(2.28) 


(2.29) 
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a 


41 


yg  {N^-mX^)  -Ni,  (  Y^-m) 
F4 


(2.30) 


a 


43 


F4 


(2.31) 


(N^-mX^)  (Y^-mu)  ~(N^-mXsU)  (Y^-m) 

F4 

(2.32) 

Y^-m 
=  ~ — 

(2.33) 

^ 

F4 

(2.34) 

Yf-mX^ 

(Y^-rn)  {N^-1^)  -(Yf-mX^)  (N^-mX^) 

(2.35) 

N,-I^ 


■^32 


iY^-m)  -  {Y^-mX^)  (N-v-mX^) 


(2.36) 


b 


41 


Yy-m 

F4 


(2.37) 
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N^-mX^ 


(2.38) 


F,=  iN^-mX^)  (  Y^-mX^)  - 


9-52  / 

aj3 

=  1. 

3-66  » 

^77 

=  -l, 

Tr 

All  other  terms  in  the  matrices  are 


(Y^-rn)  (2.39) 


zero. 
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Fig^jre  2.4  Analog  Diagram  Representing  Rudder  Dynairiics 

*  Note: 

Maximum  rudder  error  signal  ±  0.29  rad{17  degree). 
K=  0.2/sec. 

Maximum  rudder  rate  0.0698  rad/sec(3.5  deg/sec). 


C.  EXTERNAL  FORCES  AND  MOMENTS 


External  disturbances  that  act  on  a  ship  during  underway 
replenishment  operations  are  due  to  several  sources.  First, 
the  first  order  wave  forces  which  are  rapidly  varying  forces 
that  depend  on  the  seakeeping  qualities  of  the  ship  and  the 
particular  seaway  during  operations.  Slowly  varying  wave 
second  order  drift  forces  are  important  for  long  term  moving 
and  positioning  operations  while  the  same  is  true  for  current 
loads.  On  the  other  hand,  of  particular  significance  to  the 
UNREP  problem  are  the  interaction  forcesa  and  moments  caused 
by  the  proximity  of  the  two  ships.  For  this  reason  we 
concentrate  our  efforts  in  the  analysis  on  the  control  system 
performance  of  the  two  categories  of  forces:  First  order  wave 
and  interaction. 

1.  Interaction  Force  and  Moment 

These  forces  arise  mainly  due  to  changes  in  the  flow 
field  between  two  ships  in  proximity  to  one  another  and  are 
modeled  according  to  the  data  presented  in  [Ref .3] [Ref .8] .  As 
shown  in  Figure  2.6  and  Figure  2.7  these  interaction  forces 
are  functions  of  both  the  longitudinal  and  lateral  separation 
distances  of  the  two  ships.  Furthermore,  they  are  not  of 
consistent  signs  which  makes  their  effect  on  path  keeping 
control  especially  troublesome.  For  simulation  purposes  the 
data  of  these  two  figures  are  used  with  bilinear  interpolation 
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for  the  actual  values  of  longitudinal  and  lateral  ship  to  ship 
separation  distances. 

2.  First-Order  Irregular  Wave  Effects 

The  nondimens ional  first-order  irregular  sway  force 
and  yaw  moment  are  functions  of  the  regular  wave  encounter 
frequency.  The  Pierson-Moskowitz  energy  density  spectrum  is 
used  here  as  a  computer  input  to  characterize  the  sea  state. 
The  nondimens ional  force  and  moment  can  be  expressed  as [Ref .8] 


{  t)  =JcOS  ((A)^t-Cf(0)J  l-Ky(We) 


(2.40) 


WM  t)  =j*C0S  (0)gt-C„(0)g)  )^2S_^(G)g)  (2.41) 


where 


y'  (t) 
N'  (t) 

S,(wJ 
£f  (wj 

|H,(a;,)  I,  |H„(a;J  | 


is  first-order  irregular  sway  force, 
is  first-order  irregular  yaw  moment, 
frequency  of  encounter. 

Pierson  Moskowitz  energy  density  spectrum, 
phase  angle  for  sway  force, 
phase  angle  for  yaw  moment, 
the  linear  transfer  functions  at  the 
frequency  of  encounter  for  the  sway  force 
excitation  and  yaw  moment  excitation. 
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I  ^,1  ^  are  the  response  amplitude  operators  (RAO) 

for  the  sway  force  and  moment  at  the 
frequency  of  encounter. 

The  data  from  Table  2.4 [Ref. 8]  were  used  to  determine 
the  response  amplitude  operators  and  were  nondimensionlized 
with  respect  to  the  series  60  model.  For  use  with  the 
maneuvering  model  of  this  thesis  a  futher  nondimensionlization 
step  is  required. 

For  the  nondimens ional  sway  force: 

=F{(^J  x—^ -  (2.42) 


For  the  ntiidimensional  yaw  moment: 

M  ql 

H„((d^)  =7(0)^)  ^ 

•ipLiu^ 


(2.43) 


where  L,  is  the  mariner- type  ship  length. 

is  the  series  60  model  length. 

Mn,  is  the  displacement  of  series  60  model. 

The  frequency  of  encounter  is  function  of  the  absolute 
wave  frequency  u,  the  ship  speed,  the  ship- to- wave  angle  6, 
and  the  heading  angle  see  Figure  2.8.  The  equation  is 
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(2.42) 


The  Pierson-Moskowitz  Energy  Density  Spectrum  S^(w)  is 
used  to  characterize  the  sea  state.  This  is  calculated  for  a 
given  significant  wave  height  and  for  a  fully  developed,  wind- 
driven  sea.  The  characteristics  of  the  wind  speed  and  the 
corresponding  significant  wave  height  are  shown  in  Table  2.2. 
A  model  for  Pierson-Moskowitz  Spectrum  is 


where 

a=  8.1  *  lO-* 

0^  -0.74 

g=  32.174  ft/sec* 

0)=  frequency  of  the  wave  in  rad/sec 
V=  wind  speed  in  ft/sec 

The  corresponding  Pierson-Moskowitz  Spectrum  at  the 
frequency  of  encounter  is 

H-^ucosiQ-^))-^  (2.44) 
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Table  2.2  Characteristics  of  Series  60  Model 


Length  between  perpendicular 

5.0 

ft  ( 

1.32 

m) 

Beam 

0.667 

ft  ( 

0.2 

m.) 

Draft (even  keel) 

0.267 

ft  ( 

0.08 

m) 

Displacement 

33.27 

kg) 

lb(147.99 

Block  coefficient 

0.6 

Table  2.3  Characteristics  of  Wind-driven  Sea  States 


Wind  Speed 
(knots) 

' 

Sea 

State 

Significant 
Wave  Height 

Average 

Period (sec) 

10 

2 

0.6  m 

2.7 

20 

4 

2.2  m 

5.3 

30 

6 

5.0m 

8 . 0 

40 

7 

8.9  m 

10 . 7 

50 

8 

13.9  m 

13.4 

2' 


Table  2.4  Nondimensional  First-order  Sway  Force  and  Yaw  Moment  in 
Regular  Wave 


w, 

(rad/sec) 

F(aj,)  = 
sway 
force 
/(Mgr/L) 

Sway 

force 

phase 

angle 

Y  (wj  =yaw 
moment/ 

(Mgf) 

Yaw 

moment 

phase 

angle 

.293 

1.3191 

274 

. 07246 

169 

.361 

1.9473 

277 

.03741 

138 

.433 

2 . 5082 

277 

. 08306 

32 

.509 

2.9529 

280 

.22726 

19 

.588 

3.2653 

287 

.42925 

22 

.  679 

3 . 3462 

290 

.  6717 

23 

.756 

3 . 0542 

293 

. 89061 

23 

.  645 

2.364 

294 

.99661 

26 

.938 

1.5397 

293 

1.0366 

27 

1 . 034 

.6624 

281 

.9459 

24 

1.13  3 

.21465 

163 

. 73952 

19 

1.236 

.  73203 

122 

.47369 

10 

1 . 342 

1 . 033 

110 

. 19894 

399 

1.452 

1 . 0222 

101 

.  17588 

242 

1.565 

.  68012 

87 

.31972 

211 

M=  Displaced  mass 
g=  Acceleration  of  gravity 
Wave  amlpitude 

L=  Length  detween  perpendiculars 
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Table  2.5  Nondimensional  Hydrodynamic  Coefficients 


Nondimensional 

Coefficients 

Nondimensional 

Factors 

Nondimensional 
Values  (10^) 

Yv 

pL^u 

-1243 

Y^-m 

pL^ 

-1500 

Y,-m 

pL\. 

-510 

Y^-itdCg 

pL^ 

-27 

Ya 

pLV 

-270 

N. 

pL^U 

-351 

N,-mxG 

pL'‘u 

-227 

N,-I, 

pL' 

-68 

plV 

-126 

Nv 

pL'*u 

-19.7 

*  Note  : 


1.  The  primes  '  are  neglected  here  but  all 
coefficients  are  dimensionless  values. 

2.  Xg=  0,  L=length  between  perpendiculars. 

3.  u  is  the  ship  speed  15  knots, 
nondimensional  value  is  1. 

4 .  p  is  the  density  of  the  sea  water. 
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III.  CONTROL  DESIGN 


A.  CONTROL  METHODOLOGY 

The  method  for  conducting  the  Underway  Replenishemnt  at 
Sea  maneuver  requires  the  supply  ship  to  maintain  a  fixed 
speed  and  heading.  The  conning  officer  of  the  receiving  ship 
drives  it  alongside,  and  the  measurement  is  the  visual 
estimation  of  the  relative  distance  from  the  seaman's  eye. 
There  are  some  basic  measurement  methods  and  sensor  equipments 
which  are  utilized  on  the  ship  and  can  give  very  accurate 
distance  measurement.  The  information  can  be  fed  into  the 
computer  which  processes  the  data  and  automatically  maneuvers 
the  receiver  ship  to  the  desired  position. 

A  fixed  straight  line  and  a  zero  point  were  used  as  the 
course  and  the  center  of  gravity  of  the  leading  ship,  and  they 
are  the  references  for  driving  the  tracking  ship  to  a  position 
of  100  ft  from  this  line,  and  from  -524  ft  to  524  ft  with 
respect  to  this  point.  The  controller  is  at  the  tracking  ship 
(receiving  ship) .  When  the  tracking  ship  approaches  the 
leading  ship  (supply  ship) ,  the  controller  of  the  tracking 
ship  will  be  aware  of  the  excitation  due  to  the  interaction 
forces  and  moments  from  the  leading  ship.  The  dynamics  of  the 
tracking  ship  will  converge  to  the  sliding  surface  which  is 
chosen  in  order  to  drive  the  ship  to  the  desired  position. 
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At  the  station  keeping  process,  after  the  tracking  ship 
and  leading  ship  are  alongside,  a  constant  speed  which  is  a 
normalized  speed  is  considered  for  both  ships.  Also  constant 
interaction  forces  and  moments  are  exerted  at  this  stage, 
which  will  cause  a  steady  state  error.  Feedforward  control  is 
used  for  eliminating  this  error. 

The  feedforward  control  is  mainly  used  for  compensating 
the  error  of  the  ship  lateral  distance  at  steady  state.  The 
feedback  control  law  is  based  on  stability  requirements,  while 
the  feedforward  part  is  based  on  the  desired  steady  state 
accuracy.  Since  the  feedforward  control  is  a  function  of  the 
interaction  force  and  moment,  when  the  latter  are  zero,  the 
former  is  also  zero.  Therefore,  in  general  we  can  write 

Control  input  ib^)  =  Feedback  Control  {bf^) 

+Feedforward  Control  {bf^) 

During  the  control  design  stage  the  first-order  wave  force 
and  moment  and  measurement  noise  are  not  considered,  and 
perfect  and  complete  state  measurements  are  assumed.  Then  a 
kalman  filter  is  designed  for  observing  the  unmeasurable 
states  and  the  noise  filter. 

Figure  3.1  presents  a  block  diagram  of  the  optimal  control 
design  for  the  ship  replenishment  maneuver. 
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Figure  3.1  Tne  Sliding  Mode  Control  Design  for  Ship 


Replenishinent  at  Sea 


B.  SLIDING  SURFACE  DESIGN  FOR  THE  FEEDBACK  CONTROL 
1.  Reduction  of  Order 

In  sliding  mode  control,  the  equivalent  system  must 
satisfy  not  only  the  n-dimensional  state  equations  but  also  m- 
dimensional  algebraic  equations.  This  indicates  that  there 
exist  only  (n-m)  indepedent  equations  on  the  sliding  surface 
and  m  poles  are  located  at  the  origin,  where  n  is  the  number 
of  states  and  m  the  number  of  control  parameters  (1  in  our 
case)  [Ref .10]  [Ref. 12] . 

Consider  the  state  equation 

x=Ax+Bu  (3.1) 

where  states  and  matrices  refer  to  the  mathematical  model 
described  in  Chapter  II. 

The  first  state  equation  of  (3.1)  describes  the  rudder 
dynamics ; 

6  =  +  (3.2) 

and  the  rest  describe  ,  v,  r,  y,  Y,  and  N. 

The  system  of  equations  (3.1),  (3.2)  can  be  rewritten  as: 
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(3.3) 


u 


In  this  form 

X,  =  [6]  . 

V 


Y 

N 


A,,  =  [1/tj  1  by  1  matrix. 

A|2  =  0  1  by  6  matrix. 

A;,  =  Equivalent  control  6  by  1  matrix. 

Ai;  =  Equivalent  state  6  by  6  matrix. 

The  linear  switching  surface  parameters  are 
K  =  [K,  Kj] 

The  switching  surface  can  be  written  as 


p  ( x)  =K^x^  *^2^2  (3.5) 

where  K,  =1  and  is  a  1  by  6  matrix  which  is  computed  by  the 
LQR  method  for  the  equivalent  system  as  described  later. 

The  equivalent  control  is  based  on  the  system 
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x.jt)=  [1-^21  '^^21^22^2  (  ^'> 


(3.6) 


x^~-K2^K^X^  (3.7) 

with  X2  viewed  as  the  state  and  x,  as  the  control  input. 

By  substituting  (3.6)  into  (3.7),  the  reduced  order 
dynamics  can  be  rewritten  into  the  form: 

(3.8) 

Having  found  -Kj'K,  the  switching  surface  and  control  law  is 
completely  defined. 

2 .  LQR  Method  for  Determination  of  the  Sliding  Surface 

The  control  law  design  is  based  on  the  following:  A  5 
degree  of  rudder  for  path  control  is  used  when  the  heading 
deviates  5  degrees  from  zero  or  the  ship  reaches  a  lateral 
offset  of  10.43  m  (one  quarter  of  beam).  From  the  linear 
quadratic  regulator  theory  the  control  should  minimize  the 
quadratic  performance  function: 

j=lj  ix'^Ox^u^Ru)  dt  (3.9) 


where  Q  is  a  positive  semidefinite  matrix  weighting  on  the 
states,  and  R  is  a  positive  definite  control  weight. 
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The  two  weights  Q  and  R  can  be  obtained  from  the  criteria 
mentioned  above  and  can  be  evaluated  as  follows: 

Weighting  on  heading  angle  \p: 

<?22=(  5:^3)  ■'=131. 332  (3.10) 

Weighting  on  side  distance  y: 

g  -2=772. 463  (3.11) 

”  290 

Weighting  on  rudder  angle  6: 

gii=(— ^)'2  =  131.332  (3.12) 

Weighting  on  commanded  rudder  6^: 

-rii=(-^)-2=i31.332  (3.13) 

57 . 3 

The  equivalent  system  is: 

-A22^2  ^-^21^1  (3.14) 

The  Riccati  equation  is: 

A22K2-^K2A22-K2A2^Q~^aZ.K^Q=0  (3.15) 
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The  equivalent  control  is  then  given  by: 

x^  =  -Q'^aZ.K2X2  (3.16) 

3 .  System  Dynamics  and  Switches 
a.  Theory  of  Switching  Control 

Switching  control  has  been  successfully  applied  the 
nonlinear  processes  such  as  the  "Bang-Bang"  control  and  is 
widely  used  in  optimal  control.  Systems  with  switched  control 
law  represent  differential  equations  with  discontinuous  right 
hand  side.  A  problem  arises  in  that  the  conventional 
existence-uniqueness  theory  for  ordinary  differential 
equations  is  not  valid.  A  discontinuous  differential  equation 
is  defined  on  R"  as  follows:  [Ref.  10]  [Ref  .11] 

Let  S  be  defined  as  {x:p(x)=0  },  where  p  is  a  function 
from  R"  to  R,  and  S  is  (n-1)  dimensional  which  represents  the 
switching  boundary.  The  switching  dynamics  are  then  defined: 

(x=f*ix):x=pix)^0  (3.17) 

\x=f'  (x)  :  x=p  (x)  ^0 

where  and  f'  are  smooth  functions  from  R"  to  R  and  are  not 
matched  on  S  so  that  the  dynamics  are  discontinuous  at  S.  The 
mechanism  of  the  switches  are  operated  as  p  (x)  changes  sign 
which  implies  that  the  state  trajectory  passes  through  the 
surface  and  across  it. 


From  the  ideal  switching  law,  existence  of  a  sliding 
mode  requires  infinitely  fast  switching.  But  in  actual  systems 
imperfections  exist  in  the  facilities  responsible  to  the 
switching  control  such  as  delay,  hysteresis,  saturation,  etc., 
which  force  the  switches  with  a  finite  frequency.  The  ideal 
switching  can  still  be  considered  as  long  as  the  frequency  of 
the  switching  is  much  higher  than  the  dynamic  response  of  the 
system.  This  regularization  of  the  switching  dynamics  is 
schematically  depicted  Figure  3.2  and  Figure  3.3. 

The  value  p  represents  the  switching  variable:  when 
p2+l,  the  dynamics  are  described  by  +r,  and  when  ps-1,  they 
are  described  by  -r,  where  r  is  chosen  as  1 .  A  linear  relation 
of  X  and  p (x)  within  this  saturation  region  is  described  by 
the  switching  function. 

Existence  of  the  sliding  mode  requires  stability  of 
the  state  trajectory  to  the  sliding  surface  p(x)=0  and 
asymptotic  stablility  within  the  region  of  attraction  which  is 
the  region  between  +A  and  -A. 
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boundary  Layer 
Regulation 


b.  The  Lyapunov  Function 

The  second  method  of  Lyapunov  provides  a 
characteristic  function  for  analyzing  the  existence  of  the 
sliding  mode.  The  Lyapunov  function  vit,x,p)  is  continuously 
differentiable  with  respect  to  all  of  its  variabless  and 
satisfies  the  following  conditions:  [Ref. 12] 


•  v(t,x,p)  is  positive  definite  with  respect  to  p. 

•  The  total  time  derivative  of  v(t,x,p)  for  the  system  has 
a  negative  supremum  for  all  xeP  except  for  x  on  the 
switching  surface  where  the  control  inputs  are  undefined. 


The  structure  of  the  function  v(t,x,p)  is 
determined  by  the  ease  which  one  can  compute  the  switching 
parameters,  such  as  pole  placement,  Root  locus,  LQR,  or  other 
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state  space  control  theory,  for  application  of  variable 
structure  control  design.  For  all  single  input  system  a 
suitable  Lyapunov  function  is  v(t,x)  =  0.5p^(x),  which  is  a 
function  of  p  depending  on  the  control  and  has  a  negative 
derivative  function  with  respect  to  time 


1  dp‘ 

2  dt 


(3.18) 


In  the  domain  of  attraction  the  state  trajectory  converges  to 
the  surface  and  is  restricted  to  the  surface  for  all 
subsequent  times.  The  feedback  gains  associated  with  the 
optimal  design  are  computed  from  Linear  Quadratic  Regulator 
theory.  The  weighting  matrices  for  control  and  state 
associated  with  the  cost  function  are  obtained  from  (3.10)  to 
(3.13)  . 

4 .  Feedback  Control  Law 

The  method  of  determination  of  the  switching  surface 
p (x) =0  is  by  Lyapunov  function.  The  stability  function 
guarantees  that  a  sliding  mode  exist  at  every  tsto,  where  to  is 
the  time  when  the  state  trajectory  intercepts  the  surface. 

The  existence  of  the  sliding  mode  implies  that: 

•  dp (x) /dt  <  0  . 

•  p(x)=Kx  =0,  where  the  K  is  the  surface  parameters 
(feedback  gains),  and  x  is  the  state  vector. 
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From  the  chain  rule 


\ 


p  (x)  =-^x=Kx<0 
dx 


(3.19) 


we  substitute  the  constraint  equation 


x=Ax+Bu 


(3.20) 


and  the  definition 


fi  =  -r\^sign(p  (x)  ) 


(3.21) 


where  rj  is  the  weighting  for  the  switching  structure. 

This  form  guarantees  the  existence  of  a  sliding  mode  and  will 
operate  as  the  sign  of  the  state  trajectory  changes. 
Combining  (3.19),  (3.20)  and  (3.21) 

■^^^~x=-‘f]^slgn(p  ix)  )  (3.22) 

where  dp (x) /dx  =K  and  suostituting  (3.22),  then  the  control 
law  is  written  as 


u=-{ 


(3-23 

dx  dx  dx  I  ^  r  ) 


44 


There  are  two  parts  in  this  feedback  control.  They  are 
a  feedback  structure  term  for  the  first  part,  and  the  second 
part  is  for  the  switching  structure. 

In  general  the  steps  for  the  control  design  are: 

•  Determining  the  reduced  order  dynamical  equations 
governing  the  system  motion  on  the  switching  surface. 

•  Choosing  the  switching  surface  parameters  K  for  a  linear 
switching  surface  p(x)=0,  so  that  the  system  is  stat’.e  in 
the  sliding  surface  that  has  been  chosen. 

•  Augmenting  the  control  law. 

5.  Switching  Surface  Design 

Generally,  the  stability  of  the  closed  loop  system 
depends  on  the  design  of  the  switching  surface.  The  design  of 
the  switching  surface  depends  on  the  choice  of  the  surface 
parameters  and  the  dynamics  of  the  system  equation. 

The  switching  surfaces  are  designed  such  that  the 
state  trajectory  is  restricted  to  p(x)=0,  and  that  stability 
is  guaranteed. 

Again  from  the  reduced  order  dynamics  the  state 
equation  can  be  written  as: 


(3.24) 
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The  rudder  dynamics,  first  equation  of  (3.3),  will  be 
responsible  for  the  control  phase  of  the  second  equation  of 
(3.3),  which  has  an  order  n-1.  The  equivalent  control  gives 
the  equivalent  system  in  the  form 

•^2  ~-^22^2  ^■^21^1  ^  ^  -25) 


where 

Xj  is  the  (n- 1) th-order  equivalent  state  vector. 

Xi  is  the  1st -order  equivalent  control  input. 

A22  is  the  equivalent  state  matrix  (n-l)  by  (n-1) 
matrix. 

A2,  is  the  equivalent  control  matrix  (n-m)  by  m  matrix. 
The  gain  Kj  results  by  minimizing  the  cost  function 

•^2' [  [  ^2  02^2  '*'^1^2^ I  -  ^  (3.26) 

where 

i?2=[qn]  (3.27) 


Q^=diagi\q22  0  0  g^s  0  o|) 


(3.28) 
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The  weighting  factors  are  as  indicated  by  (3.10),  (3.11),  and 

(3.12),  (3.13). 

Since  (A2i,A22)  is  controllable,  it  is  possible  to 
effectively  use  classical  feedback  control  techniques  to 
compute  a  Kj  and  K,  such  that  (A22+A2iK,‘'K2)  has  the  desired 
characteristics.  Having  found  [K,  K2]  ,  the  linear  switching 
surface  is: 


[K^  K2]x=0.  (3.29) 

C.  SLIDING  SURFACE  DESIGN  FOR  FEEDFORWARD  CONTROL 
1.  Introduction 

There  are  two  types  of  disturbances  that  were 
considered  in  this  thesis.  These  are  time  varying  low 
frequency  disturbances  at  passing  process,  and  constant 
disturbances  at  station  keeping  process. 

A  steady  state  error  in  the  station  keeping  process 
indicates  that  feedback  control  is  not  enough.  The 
compensation  is  achieved  by  feedforward  control. 

The  design  objective  is  to  eliminate  the  steady  state 
error.  Prediction  of  disturbance  and  cancellation  process  is 
done  by  the  feedforward  control  generated  at  time  t.  The 
disturbances  were  assumed  to  be  the  same  at  one  time  step 
before  (t-At)  and  one  time  after  (t+At) 


47 


Ai  -r(  t+At)  -F(  t) 
aN=N(  t+At)  -N(  t) 


(3.30) 


And 


aY=Y(  t)  -Y(  t-At) 
aN=N(  t)  -N(  t-At) 


By  the  assumption  that  At  is  small  and  we  have  slowly 
varying  interaction  forces  and  moments  as  the  tracking  ship 
passes  the  leading  ship.  The  differences  will  approach  to  zero 


[lim  Ay=0 

JaC-O 

jlim  aN=0 
Ut-o 


(3.32) 


2 .  The  Feedforward  Control  Law 

The  steady  state  Conditions  are: 


r=0 

i|;  =  0=»r  =  0 

y=0=»ilf  =  -r  (3.33) 

6=0-5=6^ 
v=0 
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The  separation  distance  is  assumed  to  be  y=0,  which  is  the 
desired  position  in  which  the  tracking  ship  will  be  located. 
We  substitute  (3.33)  into  the  state  equations  and  we  get 

a^6+jbjV=Cj  (3.34) 

a^h-*-b^v=C2  (3.35) 


where 


a.^  =  Yf^(N^-mX^) (3.36) 

jb.  =  F,.  ( N^.-mX^)  -AV  ( Y^.-m)  (3.37) 

c^  =  -(Y*f^)  iY^-m)  (N+n,J  (3.38) 

a2  =  Yf,{Ni-IJ  -Nf,{Yf-mX^)  (3.39) 

b^  =  Y,.{Nf-IJ -N^AY^-mX^)  (3.40) 

c^^  =  -(Nf-lA  (Y*f.J  *(Y^-mX.)  (W+nj)  (3.41) 
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=  (3.44) 

The  results  of  these  algebraic  equations  provide  the 
amount  of  the  control  needed. 

3.  Switching  Surface  Design  For  Feedforward  Control 

As  mentioned  above  rudder  angle,  heading  angle,  and 
lateral  separation  rate  are  depend  on  the  interaction  forces 
and  moments  at  steady  state.  If  constant  disturbances  act  on 
the  ship  then  these  three  states  will  not  go  to  zero  and  a 
steady  state  error  will  be  presented.  An  appropriate 
compensation  input  will  be  found  to  eliminate  this  error. 

The  steady  state  switching  surface  is: 

p^.{x)  (3.45) 
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4 .  Augment  of  Automatic  Control  Law 

A  subsliding  surface  was  designed  for  this  feedforward 
control  input.  (An  subscript  s  indicates  that  the  states  are 
for  feedforward  calculation) 

Substitute  the  feedforward  control  switching  surface 
(3.45)  into  the  control  law  and  augment  the  reference  input  as 

b^^-{KB)  {KB)  '^r\^sign(p^)  (3.46) 

where 

6,  is  the  feedforward  control  rudder  input, 
p,  is  the  switching  surface  for  feedforward  control 
at  steady  state  condition. 
r}  is  the  weighting  factor  for  switching  control. 


0 

N 

Y 


X3  is  the  state  vector  at  steady  state. 

The  augmented  control  law  for  traclcing  design  is  completed 
and  the  schematic  block  diagram  is  shown  in  Figure  3.4. 
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inx 


c  ^ 


*  Note  for  Fig  3.4: 


fl,  =  ( N^~mx^)  F,  -  (  Y^-m)  (3.48) 

B^={Nt-I^)Y^-(Y^-mx^)N^  (3.49) 

Y^-(Y^-m)N^  (3.50 

A^=(N^-I^)  Y(,-{Y^-mx^)  Ni^  (3.51) 


fi.n,  :  1st -order  wave  force  and  moment. 
N,Y  :  Interaction  forces  and  moments. 


53 


5.  Simulation  at  Calm  Sea 

a.  Introduction 

Fortran  Program  SHIP4QQ. FOR (see  Appendix  B)  was  used 
to  simulate  the  UNREP  in  the  horizontal  plane  using  controller 
designed  by  LQR.  Simulation  is  performed  with  a  perfect  and 
full  state  feedback  at  this  stage.  Interaction  forces  and 
moments  are  calculated  by  the  bilinear  rule.  The  surge 
velocity  is  simulated  as  15  knots.  The  time  step  is  .005 
dimensionless  seconds  which  corresponds  to  10  Hz  sampling 
rate.  Two  data  files  TABLE1.DAT  and  TABLE2.DAT  contained  the 
data  points  from  Figure  2.6  and  Figure  2.7. 

b.  Simulation  Condition  Specification 

Two  types  of  ship  maneuvering  processes  are 
considered  in  this  control  test.  These  are  the  passing 
maneu'  r  and  the  station  keeping  process.  The  passing  maneuver 
is  a  process  in  which  the  tracking  ship  is  passing  the  leading 
ship  at  constant  speed  while  keeping  a  desired  course.  The 
station  keeping  process  is  a  process  in  which  interaction 
forces  and  moments  are  constant  after  both  ships  are  driven 
alongside,  where  the  same  course  and  speed  for  the  leading  and 
tracking  ship  are  maintained  at  this  stage. 

In  this  chapter,  the  sea  state  is  considered  as  a 
calm  sea.  Perfect  measurements  from  the  sensors  are  assumed. 
The  initial  conditions  for  the  tracking  ship  are  specified  as: 

•  One  nomalized  separation  distance,  y(0)=l. 
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•  Initial  heading  angle  5  degrees,  i^(0)  =  .0d7  rad. 

•  Simulation  starts  at  -4000  ft  longitudinal  distance,  and 
ends  at  8000  ft. 

•  The  leading  ship  travels  twice  as  fast  as  the  tracking 
ship  during  the  passing  process. 

The  weighting  factor  rj  in  the  switching  structure 
of  the  feedback  control  dominates  the  sensitivity  of  the 
control  law,  the  stability  of  the  state  trajectory  in  the  zone 
of  the  attraction,  and  the  capability  of  the  design  to  remain 
stable  under  external  disturbances.  Too  small  a  value 
indicates  that  the  system  may  be  slow  or  even  unstable,  too 
high  a  value  may  cause  a  chattering  problem.  The  chattering 
problem  will  not  be  discussed  here [Ref. 14] [Ref. 12]. 

D.  SIMULATION  RESULTS  AT  CALM  SEA 
1.  Passing  Process,  ijsl,  A=0.12 

Results  for  this  case  are  presented  in  Figure  3.4.  The 
four  graphs  in  this  and  all  similar  figures  show  clockwise 
from  the  upper  left  corner  the  interaction  force  and  moment, 
the  actual  rudder  angle,  the  tracking  ship  lateral  deviation 
from  the  desired  position (  100  ft),  and  the  ship  heading  in 
radians.  All  variables  are  presented  versus  the  ship- to- ship 
longitudinal  separation  with  0  corresponding  to  the  two  ships 
alongside.  It  can  be  seen  from  the  figure  that  the  control 
objectives  are  met.  The  effect  of  the  interaction  forces  is 
evident  by  the  rudder  activity  at  about  0  ft  longitudinal 
distance . 
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2.  Passing  Process,  ij=l,  A=0.06 

The  effects  of  reducing  the  boundary  layer  thickness 
are  shown  in  Figure  3.5.  The  results  are  virtually  identical 
to  those  of  Figure  3.4  which  means  that  for  the  above  range  of 
variation,  A  has  no  significant  effect.  It  should  be 
mentioned,  though,  that  this  is  true  under  the  current 
assumptions  of  calm  sea  and  perfect  feedback,  and  it  may  no 
longer  be  the  case  under  the  general  conditions  studied  in  the 
next  chapter. 

3.  Passing  Process,  i}=8,  A=0.12 

The  effects  of  increasing  the  switching  gain  17  are 
shown  in  Figure  3.6.  In  order  to  demonstrate  this  effect,  the 
commanded  rudder  angle  is  presented  instead  of  the  interaction 
forces.  It  can  be  clearly  seen  that  increasing  rj  results  in 
more  control  chattering  and  increased  overshoot  of  the 
commanded  path.  Both  of  these  are  operationally  dangerous  and 
should  be  avoided  in  practice. 

4.  Passing  Process,  17=0.5,  A=0.12 

The  effects  of  reducing  the  value  of  the  switching 
gain  are  shown  in  Figure  3.7.  The  control  response  is  now  very 
sluggish  and  the  objectives  of  the  underway  replenishment 
maneuver  are  not  met. 
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Figure  3.5  UNREP  Passing  Process,  No  Feedforward,  Calm 
Sea  (Prefect  Feedback,  tj  =  1,  A=^0.12) 
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Figure  3.5  UNREP  Passing  Process,  No  Feedforward,  Calm 
Sea(Perfect  Feedback.,  17  =  1,  A=0.06) 
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5.  Station  Keeping  Process,  rj=l,  A=0.12,  No  Feedforward 

During  the  station  keeping  process  where  the  two  ships 
travel  alongside  for  some  time,  constant  interaction  force  and 
moment  will  develop,  see  Figure  3.8.  As  a  result  in  the 
absence  of  the  feedforward  term,  a  nonzero  steady  state  error 
will  be  present.  Although  the  feedback  control  stabilizes  the 
system,  it  cannot  ensure  the  required  steady  state  accuracy. 

6.  Passing  Process,  17=1,  Aa0.12 

The  simulations  that  follow  include  the  feedforward 
term  in  the  control  lawas  developed  previously.  The  results 
are  shown  in  Figure  3.9.  Comparing  these  to  Figure  3.4.  where 
feedforward  was  not  included  we  can  see  that  for  the  passing 
process  the  effects  of  the  feedforward  term  are  very  small 
except  for  the  increased  rudder  action  at  0  ft  longitudinal 
position,  as  expected. 

7.  Passing  Process,  17=4,  A=0.06 

As  seem  from  Figure  3.10  increasing  rjand  decreasing  A 
results  in  more  path  accuracy  and  increased  control  activity. 
The  tighter  control  law  in  this  case  helps  however  during  the 
station  keeping  process  in  the  following  two  figures. 

8.  Station  Keeping  Process,  17=1,  A=0.12 

The  effect  that  the  feedforward  term  has  on  path 
accuracy  is  evident  by  comparing  the  results  of  this 
simulation  Figure  3.11,  with  those  at  the  no  feedforward  case. 
Figure  3.8.  Although  the  transient  response  during  the 
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Figure  3.9  UNREP  Station  Keeping  Process,  No  Feedforward, 
Calm  Sea  (Perfect  Feedback,  t/  =  1,  A=0.12) 
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Figure  3.10  UNREP  Passing  Process,  Calm  Sea (Perfect 

Feedback,  Disturbance  Compensation,  17  =  1,  A=0.12) 
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Figure  3.11  UNREP  Passing  Process,  Calm  Sea (Perfect 

Feedback.,  Disturbance  Compensation,  t}  =  4,  A=0.06) 


approach  phase  is  identical,  the  steady  state  path  deviation 
becomes  much  smaller  and  asymptotically  reaches  zero. 

9.  Station  Keeping  Process,  i;=4,  A=0.06 

The  results  of  this  simulation  are  shown  in  Figure 
3.12  where  we  can  clearly  see  the  increased  rudder  action  and 
path  overshoot.  These  result  from  increasing  rj  and  decreasing 
A  as  in  the  passing  process. 
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Figure  3.12  UNREP  Station  Keeping  Process,  Calm  Sea (Perfect 
Feedback,  Disturbance  Compensation,  17  =  1,  A=0.12) 
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Figure  3.13  UNREP  Station  Keeping  Process,  Calm  Sea (Perfect 
Feedback,  Disturbance  Compensation,  17=4,  A=0.06) 


67 


E.  REMARKS 


A  sliding  mode  controller  has  been  designed  for  the  UNREP 
which  utilized  a  LQR  method.  The  controller  have  been  proved 
to  be  very  effective  for  a  time  varying  disturbance  (passing)  . 
Disturbance  compensation  was  achieved  by  means  of  a 
feedforward  term  in  the  automatic  controller. 

The  results  demonstrate  excellent  tracking  and  path 
keeping  characteristics  in  the  presence  of  strong  interaction 
forces  and  moments  without  any  compromise  in  the  stability  and 
robustness  properties  at  calm  seas. 

In  this  control  test  results  all  states  are  measurable, 
but  in  reality  for  a  ship  operating  at  sea  not  all  states  can 
be  measured.  So  a  Kalman  filter  design  is  needed  as  done  in 
Chapter  IV. 


68 


IV.  KALMAN  FILTER 


A.  INTRODUCTION 

1.  Disturbances 

As  was  demonstrated  in  the  previous  chapter  the 
optimal  sliding  mode  control  law  was  able  to  guarantee 
stability  and  steady  state  accuracy.  For  that  it  required 
knowledge  of  all  states  including  the  interaction  forces  and 
moments.  An  observer  is  therefore  necessary  in  order  to 
provide  an  estimate  as  accurate  as  possible  of  all  system 
states.  The  first  order  wave  forces  can  be  effectively  modeled 
as  white  noise  components  since  they  vary  very  rapidly 
compared  to  the  dynamics  of  the  maneuvering  ship.  A  shaping 
filter  is  introduced  in  order  to  model  the  interaction  force 
Y  and  moment  N  as  exponentially  correlated  noise  driven  by 
white  noise.  In  other  words 


N^-  —  N+~W, 


N 


‘W 


Y= 


(4.1) 


where  Wj^,  Wy  are  white  noise  components  and  Tn  and  Ty  are  time 
constants  of  the  filter  which  are  approximately  equal  to  one 
(the  tine  that  it  takes  to  travel  ine  ship  length) . 
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The  power  spectral  density  of  Y  and  N  is 


(4.2) 


where  and  Qy  are  estimated  using 


Ok  ^ 


(4.3) 


Qy=2(s\Ty  (4.4) 

The  root  mean  aquare  values  of  N  and  Y  denoted  by  and 
Cy  can  be  estimated  using  the  typical  curves  for  the 
interaction  forces  shown  in  Figure  4.1.  This  calculation  gives 


548x10’®  0 

0  8.97x10’® 


(4.5) 


2 .  Measurement  Noise 

The  measurable  states  in  this  study  are  heading 
angle  (!/•),  yaw  rate(r),  and  lateral  separation  distance  (y). 
The  output  equation  is 
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0 

1 

0 

0 

0 

0 

0 

•K 

= 

0 

0 

0 

1 

0 

0 

0 

.y 

0 

0 

0 

0 

1 

0 

0 

6 

V 

^1 

yi 

N 

Y 


(4.6) 


where  is  the  output  vector,  the  output  matrix 


0 

1 

0 

0 

0 

0 

0 

c= 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

and  V,  is  the  measurement  noise.  The  power  spectral  density  R^, 
of  Vn  can  be  evaluated  using  the  numerical  values  in  Table 
4.1.  The  weighting  matrix  R,  is  written  in  the  fomn 


rv^j  0  0 

0  rv-^  0 

0  0  rvjj 


(4 . 7) 
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Table  4.1  The  Power  Spectral  Density  of 
Noises 


Measurement 


Measurement 

Source 

a 

T 

gyro 

compass 

rvj,  = 

4.612*10-® 

r  .2 

0.1 

sec 

r 

rate  gyro 

rv22  = 

3.218*10'^ 

0 . 1 
/s 

0.1 

sec 

V 

radar 

rv33  = 

4.502*10'* 

3.3 

ft 

0.1 

sec 

Note: 


o:  Random  standard  deviation  fr'-r  the  true 
value (  units  are  degrees  for  Heading  angle 
and  yaw  rate,  ft  for  y) . 
t:  The  correlation  time  . 

All  values  of  are  :jn  rondimensional 
values . 
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B.  DESIGN  OF  KT^LMAN  FILTER 


In  order  to  avoid  enlarging  the  plant  noises  and 
measuremant  noises,  an  optimal  stochastic  observer  is 
required.  The  state  vector  and  the  covariance  error  are 
described  as  follows 


£’[^(  to)  ]  =Xo 

£■[  [x(  to)  -Xol  [x(  to)  -Xq]  T  =Po 


(4.8) 


The  estimate  of  the  state  at  to  is  assumed  to  be 

x(to)=3co  (4.9) 


The  disturbance  in  the  state  equations  has 

ElwiO]  *0 

E[w(  t)  w"(  t)  j  =0(t)  b  {t-z) 


(4.10) 


The  measurement  noise  has 

E[  v(  t)  ]  =0 

E[v(t)  v~(  t)  j  =E(  t)  6  ( t-T) 


(4.11) 


Design  of  the  Kalman  Filter  is  based  on  the  cost  function 


(4.12) 


J=^  [  (a'q-Xq)  ^Pq*^  (Xq-Xq)  ] 

+  ^  f  "  [w'^0'^w+  (y^-Cx)  (y^-Cx)  ]  dt 

2  J  to 

where  the  first  term  minimizes  the  error  in  the  initial 
estimate,  the  second  term  minimizes  the  effect  of  w,  and  the 
third  term  minimizes  the  error  in  the  estimate  of  x. 

The  Kalman  Filter  has  the  familiar  form 

^=Ax^B^u+L(y^-Cx)  (4.13) 


L  is  the  time-varing  Kalman  gain  matrix, 

L=PC^R-'^ 


(4.14) 


P  is  the  solution  of  the  Riccati  differential  equation, 
P=AP^PA  '^^B^QB2-PC  ’^R-'^CP  (4.15) 

Bj  is  the  disturbance  matrix;  and 

P(to)  =Po. 

At  steady  state  P  is  the  solution  to  the  algebraic  Riccati 
equation 

AP+Pa’^->-B2QB2+PC'^R-'^CP=0  (4.16) 
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In  such  a  case  the  Kalman  filter  gain  matrix  is  time- 
invariant  and  the  Kalman  filter  is  a  full  order  observer  with 
its  gains  tuned  in  an  optimal  way  to  the  system  at  hand. 
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Figure  4.1  Turning  Phase  UNREP  Passing  Process,  Open  Sea, 
Disturbance  Compensation  7j  =  l,  A=0.12 
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Figure  4.2  Turning  Phase  UNREP  Passing  Process,  Open  Sea, 
Disturbance  Compensation  77=4,  A=0.06 
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C.  SIMULATION  AT  OPEN  SEA 


1.  Introduction 

In  the  following  we  present  results  for  both  the 
passing  and  the  station  keeping  process.  The  sliding  mode 
linear  quadratic  Gaussian  compensator  is  used  along  with  the 
feedforward  term  developed  in  chapter  III.  All  results  are 
based  on  UNREP  in  a  random  sea  of  significant  wave  height  of 
16  ft.  Solid  curves  in  the  graphs  denote  the  actual  values  of 
the  variables  while  dotted  curves  correspond  to  the  estimated 
values  from  the  Kalman  filter. 

2.  Passing  Process, A=0. 12,  t}=1 

The  results  of  the  simulation  are  shown  in  Figure  4.1. 
Convergence  to  the  desired  track  is  achieved.  The  steady  state 
deviation  of  the  tracking  ship  from  the  commanded  path  is 
minimal.  The  steady  state  use  of  the  rudder  angle  is  also  very 
satisfactory. 

3.  Passing  Process,  j;=4,  A=0.06 

The  results  of  the  simulation  are  shown  in  Figure  4.2. 
Initial  conditions  were  selected  on  the  desired  path  so  this 
simulation  demonstrates  the  effects  of  the  passing  ship 
disturbance.  The  two  ships  are  located  beam-to-beam  at  0  ft 
longitudinal  distance.  The  maximum  deviation  of  the  tracking 
ship  is  ±  3  ft. 
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Figure  4.3  Turning  Phase  Station  Keeping  Process,  Open  Sea, 
Disturbance  Compensation  17=1,  A=0.12 
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4.  Station  Keeping  Process,  ti=l,  AsO.12 

The  results  of  this  simulation  are  shown  in  Figure 

4.3.  Steady  state  accuracy  is  maintained.  There  appears  to 
exist  an  error  in  the  estimates  of  Y  and  N  despite  the  zero 
steady  state  error  in  the  path  deviation.  This  is  attributed 
to  the  existence  of  the  first  order  wave  force  and  moment. 
Since  there  is  no  explicit  estimation  scheme  for  these,  their 
effect  is  effectively  combined  by  the  observer  into  an 
equivalent  interaction  force  and  moment  which  of  course  differ 
from  the  actual. 

5.  Station  Keeping  process,  i}=4,  A=0.06 

The  results  of  this  simulation  are  shown  in  Figure 

4.4.  Compared  to  the  simulation  shown  in  Figure  4.3  we  can  see 
the  effects  of  the  tighter  control  law  in  this  case.  Both  the 
transient  and  steady  state  path  deviation  are  reduced  while 
the  rudder  activity  is  slightly  increased. 

D.  CONCLUSIONS  AND  REMARKS 

The  sliding  mode  Linear  Quadratic  Gaussian  compensator  was 
able  to  adequately  control  the  lateral  motion  of  a  surface 
ship  during  underway  replenishment  in  an  open  sea.  The 
feedforward  term  ensured  state  path  accuracy  in  the  presence 
of  interaction  forces  and  moments  as  well  as  first  order  wave 
effects.  The  Kalman  filter  was  able  to  minimize  the  effects  of 
the  measurement  noises  and  external  random  disturbances  while 
at  the  same  time  it  provided  the  controller  with  an  accurate 
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estimate  of  the  states  needed  for  control.  Different  control 
options  were  considered  with  regards  to  their  transient  and 
steady  state  characteristics.  These  comparisons  provided 
insight  towards  appropriate  selection  of  the  boundary  layer 
thickness  and  switching  gain  of  the  sliding  mode  control. 

Suggestions  for  further  research  include  the  application 
of  integral  control  instead  of  disturbance  estimation  and 
compensation  for  eliminating  the  steady  state  path  error. 
Furthermore,  the  effects  of  constant  disturbances,  such  as 
currents,  and  slowly  varying  second  order  wave  draft  forces 
should  be  investigated. 
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APPENDIX  A.  SUMMARY  FOR  COMPUTER  PROGRAMS 


The  complete  listing  of  the  computer  programs  written  and 
used  in  this  thesis  are  in  APPENDIX  B  and  APPENDIX  C. 

The  MATLAB  program  in  APPENDIX  B  is  used  for  calculating 
feedback  gains  and  observer  gains  using  the  LQR  method,  which 
are  utilized  for  sliding  surface  design  and  Kalman  filter 
design. 

The  FORTRAN  program  in  APPENDIX  C  is  a  complete  listing 
for  PASSING  PROCESS  simulation.  The  modified  program  for  the 
STATION  KEEPING  PROCESS  is  not  listed  here,  the  only  change  is 
the  following  FORTRAN  code  logic  expression: 

if  (b.ge.O)  b=0. 
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APPENDIX  B.  MATLAB  LQR  DESIGN  PROGRAM  FOR  UNREP  SLIDING  MODE 
CONTROLLER  AND  KALMAN  GAINS 


mass=0 . 0 ; 
xg=0. 0; 
iz=0 . 0 ; 
tr= . 2 ; 

nvdot=- 19 . 7e- 5 ; 
xudot=- SoOe- 5 ; 
xu=- 120e- 5 ; 
yv=-1243e-5; 
yvdot= - ISOOe- 5 ; 
yr=- 510e- 5 ; 
yrdot=-27e-5; 
ydr=270e- 5 ; 
nv= - 351e- 5 ; 
nr=-227e-5; 
nrdot= - 68e- 5 ; 
ndr=- 126e- 5 ; 
xn=4 . 62e- 5 ; 
yn= - . 52e- 5 ; 
nn= . 26e- 5 ; 

apl =yvdot  *nrdot - vrdot  *nvdot 
a32=  (yrdot*r:-'^  y  v*nrdot ) /apl  ; 
a33=  (yrdot'*nr-yr*nrdot )  /apl; 
a35=yrdot/apl ; 
a36=^-nrdot/apl  ; 
ap2=nvdot*yrdot -nrdot*yvdot 
a42= (nv*yvdot - nvdot*yv) /ap2 
a43= (nr*yvdot-yr*nvdot ) /ap2 
a45=yvdot/ap2 
a4  6  =  - nvdot /ap2 

b21= (ndr*yrdot-ydr*nrdot) /apl ; 
b3i= (ndr*yvdot-ydr*nvdot) /ap2 

a=[0  0  1  0  0  0;0  a32  a33  0  a35  a36;0  a42  a43  0  a45  a46;  ... 

110  000;0000-10;00000-l]; 

b= [0  b21  b31  0  0  0] ' ; 
aab=b31  'a32-b21*a42 ; 
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aal= (a42*a35-a45*a32) /aab; 

aa2 =a3  5  *a42 /aab ; 

bbl= (a35*b31-a45*b21) /aab; 

bb2  =a3  6  *a42 /aab ; 

tau= . 2 ; 

ad=l  ; 

all= [ - 1/ tau] ; 

al2=  [0  0  0  0  0  0]; 

bb= [1/tau  000000]'; 

aa=  [all  al2 ;b  a] ; 

q=diag(  [131.3  0  0  771.69  0  0]); 

qq=diag(  [131.3  131.3  0  0  771.69]); 

r= [131.3]  ; 

c=[l  0  0  0  0  0;0  0  1  0  0  0;0  0  0  1  0  0]; 

[k,s] =lqr (a,b,q,r) ; 
ss= [1  k] 

gg=inv (ss*bb) *ss*aa 
kn=inv {ss*bb) ; 
ca=aa-bb* (gg-knl) ; 

CC=[0  100000  ;  0001000;  000010  0]; 

gl=[0  000010  ;  000000  1]'; 

ql=[1.548e-8  0;0  8.97e-8]; 

rr=diag ( [4.612e-8  3.213e-7  4.502e-6]); 

l=lqe (aa, gl , cc , ql , rr) 

f =aa- l*cc*aa 

bf =bb- l*cc*bb 
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APPENDIX  C.  UNREP  SIMULATION  PROGRAM:  SLIDING  MODE 
CONTROLLER  DESIGNED  BY  LQR  METHOD 


program  shipdynamics 

c  Simulation  program  for  UNREP 
c 

c  Written  by  LT  FU,  HSU-SHENG  SIMON 
c 

c  The  program  controls  one  rudder  input  with  a  sliding  mode 
c  design  for  path  keeping  and  tracking, 
c 

c  Speed  is  constant  l  normalized  value, 
c 

c  The  initial  conditions  for  ship  are  follows:  1 
c  nondimensional  leteral  seperation  distance,  -4000  ft 
c  longitudinal  distance,  0.087  rad  heading  angle, 
c 

c  The  disturbances  include  interaction  forces  and  moments, 
c  and  first-order  wave  sway  force  and  yaw  moment, 
c 

c  The  observed  states:  yaw  rate,  lateral  separation  rate, 
c  interaction  forces  and  moments, 
c 

c  The  measurements:  heading  angle,  lateral  separation 

c  distance,  and  rudder  angle, 

c 

c  The  outputs  of  the  program  are  written  to  files,  these 

c  files  are  then  accessed  by  the  MATLAB  to  generat  output 

c  graphics, 

c 

c  This  program  utilized  subroutines  for  different  functions 
c  that  are  required: 

c  subroutine  locate:  find  the  interaction  forces  and 

c  moments  at  the  real  time  and  use 

c  subroutine  trlint  as  interpolation, 

c  subroutine  trlint:  incooprate  with  subroutine  locate. 
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c  subroutine  randm  :  generate  random  signal, 

c  subroutine  trap  :  Trapzoid  integration, 

c  subroutine  first  :  first  order  wave  calculation, 

c 

c  Two  data  files  for  interaction  forces  and  moments: 
c  Tablel.dat:  interaction  forces, 

c  Table2.dat:  interaction  moments. 

real  mass , nrdot , nvdot , nr , nv, ndr , nn 
real  kl , k2 , k3 , k4 , k5 , k6 , k7 , kn , i z , mo , a , b , 1 
real  w, coef , fl,nl,wn,wy,we2 (5) ,ws 

real  ol (7 , 3 ) , ae (7 , 7) ,ba (7, l) , be (7 , 2 ) ,bd (7, 2) ,sf2 (5) , 

1  ym2  (5) 

real  time 
c 

dimension  aa (23 ) , bb (11) , tablel (23,11) , table2 (23,11) 
data  aa/-550, -500, -450, -400, -350,-300, -250, -200, -150, 
1  -100, -50,0,50,100,150,200,250.300,350,400, 

1  450,500,550/ 

data  bb/50, 60, 70, 80,90,100, 110, 120, 130, 140, 150/ 
data  01/0,2.6383, -1.4846,6.6537, .165, .0167, - .0133, 

1  0, .9536, -3.8899,21.5836, - .911, .1926, - .1238, 

1  0, .0017,3.3203, - .0651,2.566, .0107, .0507/ 

C 

m  =23 

n  =11 

1  =527.8 

iseed  =17 


c 


mass  = 

xg 

iz 

nvdot =- 


0.0 


0.0 


0.0 

19.7*l.e-5 


tr  =  .2 

xudot=-  850,0*l.e-5 
XU  =-  120.0*l.e-5 


yv  =-1243 .0*1. e-5 
yvdot=- 1500 . 0*1 .e-5 
yr  =-  510. 0*1. e-5 
yrdot=-  27. 0*1. e-5 
ydr  =  270.0*1.6-5 
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c 

c 

c 


70 

60 


nv 

nr 

nrdot=- 


ndr 

xn 

yn 

nn 


351.0*1.6-5 

227.0*1.6-5 

68.0*1.6-5 

126.0*1.6-5 

4.62*1.6-5 

.52*1.6-5 

.26*1.6-5 


d6fine  syst6m  matrix 


do  60  i  =1,7 
do  70  j  =1,7 
a6 (i, j ) =0 . 0 
ba (i , 1) =0 . 0 
be (i, 1) =0 . 0 
bd{i,l)=0.0 
be ( i , 2 ) =0 . 0 
bd(i,2) =0.0 
eontinu6 
eontinu6 

ael  =  nrdot*yvdot-yrdot*nvdot 

a62  =  nvdot*yrdot-nrdot*yvdot 

ae (1, 1) =-l/tr 
a6(2,4)=  1. 

a6 (3 , 1) =  (ndr*yrdot-ydr*nrdot) /aal 

a6(3,3)=  (yrdot*nv-yv*nrdot) /aal 

a6(3,4)=  (yrdot*nr-yr*nrdot) /aal 

a6(3,6)=  yrdot/ael 

aa (3 , 7) =-nrdot/ael 

aa (4 , 1) =  (nar*yvdot-ydr*nvdot) /a62 

ae(4,3)=  (yvdot*nv-yv*nvdot ) /a62 

a6(4,4)=  {nr*yvdot-nvdot*yr) /a62 

a6(4,6)=  yvdot/ae2 

aa (4 , 7) =-nvdot/a62 

aa (5,2) =  1. 

a6(5,3)=  1. 

aa (6, 6) =-l. 

a6{7,7)=-l. 

ba(l,l)=  1/tr 

be (6 , 1) =  1 . 

be (7,2)=  1. 
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bd(3,l)=  yrdot/ael 
bd (3 , 2) =-nrdot/ael 
bd(4,l)=  yvdot/ae2 
bd (4 , 2) =-nvdot/ae2 
c 

open (10 , f ile= ' sim.dat ' , status=' old' ) 
open (11 , f ile= ' siml . res' , status= ' new' ) 
open (12 , f ile= ' sim2 . res ' , status= ' new' ) 
open ( 13 , f ile= ' sim3 . res ' , status= ' new' ) 
open (14 , f ile= ' tablel.dat ' , status= ' old' ) 
open (15 , f ile=' table2 .dat ' , status= ' old' ) 
open  (16 ,  f  ile= '  siin4  .  res'  ,  status= '  new' ) 
open (18 , f ile= ' sim5 . res ' , status=' new' ) 
open ( 19 , f ile= ' sim6 . res ' , status= ' new' ) 
open (20 , f ile= ' simol . res ' , status=' new' ) 
open (21 , f ile= ' simo2 . res ' , status=' new' ) 
c 

clOO 

do  40  i=l,m 

read(14,30) (tablel (i, j ) , j=l,n) 
read (15,30) (table2 (i,j) ,j=l,n) 

30  format (Ilf 7. 2) 

40  continue 

c 

c 

readdO,  *)  stime, delta 
read (10 ,*)sl,s2,s3,s4,s5,s6,s7 
read ( 10 , * ) kl , k2 , k3 , k4 , k5 , k6 , k7 , kn 
read (10, *) sigsat,eta,etal,psi,y,xref ,yref ,ws 
c 

ws=ws* .51444/, 3048 
c 

c  define  state  v, r,dr,drc,psi,y 

c 

is im=s time/delta 

V  =0.0 

u  =1.0 
r  =0.0 
dr  =0.0 
drc  =0.0 
X  =0.0 
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y  =0.0 
yf  =0.0 
mo  =0.0 
c 

c  estimated  states 

c 

odr  =0.0 
opsi  =0.0 
ov  =0.0 
or  =0.0 
oy  =0.0 
oyf  =0.0 
omo  =0.0 
c 

c  simulation  start 

c 

do  1  i=l,isim 
if  {b. ge . (3 . 5*xref ) )  stop 
c 

c  First  order  wave 

c 

call  f irst (psi, ws, time, f 1, nl) 
c 

c  White  noise 

c 

call  randm ( iseed, wn) 
call  randm ( iseed, wy) 
wn= (wn*l . e- 5 ) 
wy= (wy*l . e- 5 ) 
c 

c  Interaction  Force  and  Moment 

c 

a  =yref+y*l 
b  =x*l-xref 
call  locate (bb, n, a, nx) 
call  loctae (aa,m,b,my) 

call  trlint {m,n,my,nx, tablel,bb,aa,a,b,yf ) 
call  trlint (m,m,my,nx, table2,bb,aa,a,b,mo) 
yf=yf*l.e-5 
mo=mo*l.e-5 

if  (b.lt . (-550) .or.b.gt.550)  then 
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yf  =  0.0 
mo=0 . 0 
endif 
c 

al =mas s - yvdot 

bl=mass*xg-yrdot 

cl=yv*v+ (yr-mass) *r+ydr*dr+yf +f 1 

a2  =ina  s  s  *  xg  -  nvdo  t 

b2=iz-nrdot 

c2=nv*v+ (nr-mass*xg) *r+ndr*dr+mo+nl 
a3=mass-xudot 
c 

drdot  =-dr/tr+drc/tr 
psidot=  r 

vdot  =  (cl*b2-c2*bl) / {al*b2-a2*bl) 
rdot  =  {cl*a2 - c2*al) / (a2*bl-al*b2 ) 
ydot  =  sin (psi) +v*cos (psi) 
xdot  =  cos (psi) -v*sin (psi) 
yfdot  =-yf+wy 
modot  =-mo+wn 
c 

c  Euler  Integration  for  Actual  States 

c 

dr  =dr  +drdot  *delta 
psi=psi+psidot*delta 
V  =v  +vdot  *delta 
r  =r  +rdot  *delta 
y  =y  +ydot  *delta 
X  =x  +xdot  ♦delta 
yf  =yf  +yfdot  *delta 
mo  =mo  +modot  *delta 
c 

c  Measurement  Noises 

c 

call  randm ( iseed, rnpsi) 
call  randm (iseed, rnr) 
call  randm (iseed, rny) 
rnpsi=rnpsi*0 . 0035 
rnr  =rnr  *0.0002 
rny  =rny  *0.0063 
psime=psi  +rnpsi 
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rme  =r  +rnr 

yme  =y  +rny 

c 

c  Feedforward  Control  Rudder  Angle  Design 

c 

aal  =  ydr*nvdot-ndr*yvdot 

bbl  =  nvdot*yv-yvdot*nv 

ccl  =- (oyf +f 1) *nvdot+ (omo+nl) *yvdot 

aa2  =  ydr*nrdot-ndr*yrdot 

bb2  =  yv*nrdot-nv*yvdot 

cc2  =  (omo+nl) *yrdot- (oyf +fl) ♦nrdot 

drcO  =  (ccl*bb2-cc2*bbl) / (aal*bb2-aa2*bbl) 

vO  =  (ccl*aa2-cc2*aal) / (bbl*aa2-bb2*aal) 

psiO  =-v0 

sigmal=sl*drc0+s2*psi0+s3*v0+s6*omo+s7*oyf 
if  (abs (sigmal) . It . sigsat )  satsgnl=sigmal/sigsat 
if  (sigmal. le. -sigsat)  satsgnl=-1.0 

if  (sigmal. gt.  sigsat)  satsgnl=  1.0 

aO  =-etal**2*kn*satsgnl 

drO  =  drc0+ (kl*drc0+k2*psi0+)c3*v0+k6*omo+k7*oyf ) -aO 
c 

c  rudder  input  calculation 

c 

sigma=sl*odr+  (s2'*opsi+s3*ov+s4*or+s5*oy+s6*omo+s7*oyf ) 
if  (abs (sigma) . It . sigsat )  satsgn=sigma/sigsat 
if  (sigma . le .- sigsat )  satsgn=-1.0 

if  (sigma. ge.  sigsat)  satsgn=  1.0 

c 

drchat=- (kl*odr+k2*opsi+k3*ov+k4*or+k5*oy+k6*omo+ 

1  k7*oyf) 

drcbar=  eta**2*kn*satsgn 
drc  =  drchat-drcbar+drO 

drf  =  drchat-drcbar 

c 

if  (drc.ge.  0.4)  drc=  0.4 
if  (drc . le . - 0 . 4 )  drc=-0.4 
c 

c  Observed  States 

c 

obl=ae (1,1) *odr+ae (1,2) *opsi i ae ( 1 , 3 ) ♦ov+ae (1,4) *or+ 

1  ae ( 1 , 5 ) *oy+ae (1,6) *omo+ae (1,7) *oyf 
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ob2=ae (2 , l) *odr+ae (2,2) *opsi+ae (2,3) *ov+ae (2,4) *or+ 
ae (2 , 5) *oy+ae (2,6) *omo+ae (2,7) *oyf 
ob3=ae (3,1) *odr+ae (3,2) *opsi+ae (3,3) *ov+ae (3,4) *or+ 
1  ae ( 3 , 5 ) *oy+ae (3,6) *omo+ae (3,7) *oyf 

ob4=ae (4,1) *odr+ae (4,2) *opsi+ae (4,3) *ov+ae (4,4) *or+ 
1  ae (4, 5) *oy+ae (4,6) *omo+ae (4,7) *oyf 

ob5=ae (5,1) *odr+ae (5,2) *opsi+ae (5,3) *ov+ae (5,4) *or+ 
1  ae (5 , 5) *oy+ae (5,6) *omo+ae (5 , 7) *oyf 

ob6=ae (6,1) *odr+ae (6,2) *opsi+ae (6,3) *ov+ae (6,4) *or+ 
1  ae (6, 5) *oy+ae (6,6) *omo+ae (6,7) *oyf 

ob7=ae (7,1) *odr+ae (7,2) *opsi+ae (7,3) *ov+ae (7,4) *or+ 
1  ae ( 7 , 5 ) *oy+ae (7,6) *omo+ae (7,7) *oyf 

ob23  =bd(3,l)*(nl)+bd(3,2)*fl 
ob24  =bd(4, 1) * (nl) +bd(4,2) *fl 

odrdot  =obl+ol (1,1)* (psime-opsi) +ol (1,2) * (rme-or) + 

1  ol (1, 3) * (yme-oy) +ba(l, 1) *drc 

opsidot=ob2+ol (2,1)* (psime-opsi) +ol (2,2)* (rme-or) + 

1  ol (2 , 3 ) * (yme-oy) 

ovdot  =ob3+ol (3,1)* (psime-opsi) +ol (3,2)* (rme-or) + 

1  ol (3 , 3) * (yme-oy) +ob23 

ordot  =ob4+ol (4,1)* (psime-opsi) +ol (4,2) * (rme-or) + 

1  ol (4, 3) * (yme-oy) +ob24 

oydot  =ob5+ol (5,1)* (psime-opsi) +ol (5,2) * (rme-or) + 

1  ol  (5 , 3 )  *  (yme-oy) 

omodot  =ob6+ol (6,1)* (psime-opsi) +ol (6,2)* (rme-or) + 

1  ol  (6, 3)  *  (yme-oy) 

oyfdot  =ob7+ol (7,1)* (psime-opsi) +ol (7, 2) * (rme-or) + 

1  ol  (7, 3)  *  (yme-oy) 

c 

c  Eular  Integration  for  Estimate  States 

c 

odr  =odr  +odrdot  *delta 
opsi=opsi+opsidot*delta 


ov 

=ov 

+ ovdot 

*delta 

or 

=or 

+ordot 

*delta 

oy 

=oy 

+oydot 

*delta 

oyf 

=oyf 

+oyfdot 

*delta 

omo 

=omo 

+omodot 

*delta 

time=i*delta 

if  (mod ( i , 20) . eq. 0)  then 
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write (11,20) time , psi , v, y , r , dr 

format (fl0.4,lx,fl0.4,lx,fl0.4,lx,fl0.4,lx,fl0.4, 

1  lx,fl0.4) 

write (12,*)x,y,a,b 
write (13 , *) drchat , drcbar, drc, drf 
write (16, *)yf ,mo,dr0 
write (19 , *) wy , wn 
write (20,50) odr , opsi, ov, or, oy 
50  format (fl0.4,lx,fl0.4,lx,fl0.4,lx,fl0.4,lx,fl0.4) 

write (21, *) omo, oyf 
write (18, *) f l,nl 
end  if 

1  continue 

stop 
end 

c303 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  locate (xx, n, x, j ) 
c 

c  This  is  a  subroutine  for  finding  the  interaction  forces 
c  and  moments.  For  location  the  interaction  forces  and 
c  moments  in  the  input  file  TABLE1.DAT  and  TABLE2.DAT 

c  correspond  the  lateral  distance  and  longitudinal 

c  distance, 

c 

dimension  xx(n) 
c 

do  20  i=l,n 
if  (x.eq.xx(i)  )  j=i 
20  continue 

do  10  i=l,n-l 

if  (x.gt .xx(i) .and.x.lt.xx(i+l) )  then 

j  =  i 

end  if 

10  continue 

return 
end 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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c 

subroutine  trlint (m, n, i , j ,ya, xla, x2a, xl , x2 , y) 
c 

c  This  subroutine  is  used  for  calculating  the  interaction 

c  forces  and  moments,  the  data  get  from  SUB  LOCATE  and  the 

c  values  calculated  by  interpolation, 
c 

dimension  ya (m, n) , xla (n) , x2a  (m) 
real  t,u,tl,ul 
c 
c 

yl=ya(i  ,j  ) 
y2=ya(i  ,j+l) 
y3=ya(i+l,j  ) 
y4=ya (i+1 , j  +1) 
c 

t= (xl-xla ( j ) ) / (xla { j  +1) -xla ( j ) ) 
u= (x2-x2a (i) ) / (x2a (i+1) -x2a (i) ) 
c 

tl=yl+u* (y3 -yl) 
ul=y3+u* (y4-y2) 
y  =tl+t* (ul-tl) 
c 

return 

end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  randm ( iseed, rndm) 
c 

c  generate  a  gaussian  random  number 
c  iseed  -  input  a  integer  variable  seed 

c  output  seed  is  changed  during  execution  for  next 

c  call 

c  rndm  -  gaussian  random  variable  with  zero  mean  and  unit 

c  variance 

c 

real  rndm, rndl , rnd2 
integer  iseed, jr,kr,mr 
parameter  ( j r=5243 , kr=55397 ,mr=262l39 ) 
iseed=mod ( iseed* j  r+kr,mr) 
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rndl  = (iseed+ . 5) /mr 
iseed=mod ( iseed* j  r+kr , mr) 
rnd2  = (iseed+ . 5) /mr 

rndm  =sqrt ( -2 . * log (rndl) ) *cos (6 . 2831853 *rnd2 ) /2 . 5 
return 
end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  trap (nf ,vx,vy, out ) 
c 

c  numerical  integration  routine  using  the  trapezoidal  rule 
c 

dimension  vx(l),vy(l) 
real  out 
c 

nl=nt- 1 
out  =  0 . 0 
do  1  i=l,nl 

outl=.5* (vx(i) +vx(i  +  l) ) * (vy  (i  +  1) -vy (i) ) 
out  =out+outl 
1  continue 

return 
end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  f irst (psi , ws , time, f 1 , nl) 
c 

c  This  subroutine  is  used  for  calculating  the  sway  force 

c  and  yaw  moment  component  from  first-order  wave 

c  excitation, 

c 

c  The  inputs  of  the  subroutine  are: 

c  (1)  The  encounter  frequency (wj  : 

c 

c  Wj=w- (w''2 ) /g*u*cos  ( theta-psi) 

c  w  is  the  actual  wave  frequency, 

c  theta (0)  is  the  wave- to- ship  angle, 

c  psi(\l/)  is  the  ship  heading, 

c  g  is  the  acceleration  of  gravity, 

c  u  is  the  ship  forward  speed(15  knots), 

c 

c  (2)  The  wind  speed (ws). 
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c 

c  (3)  The  RAO  of  sway  force  and  RAO  of  yaw  moment  are  the 

c  function  of  encounter  frequency  . 

c  (4)  The  Pierson-Moskowitz  spectrum  is  the  function  of  we 
c  and  wind  speed  . 

c  The  associated  force  and  moment  are  calculated  from 

c  Trapizoid  Intergration  and  be  the  outputs  of  the 

c  subroutine. 

c  The  phase  angles  correspond  to  the  sway  force  and 

c  moment  are  generated  by  random  process  (irregular  wave) 
c  (epsf,epsm)  from  0  to  2tt 

c397 

real  sf (19) ,ym(l9) ,wee (19) ,wf (19) , wm(l9) , vecx(19) , 

1  vecm(19) 

real  we, alpha, beta, ws , g, epsf , epsm, sx ( 19 ) , psi , wl (19 ) , 

1  w2(19) 

real  raom. (19)  ,  raof  (19)  ,  time, model,  1,  Im,  rho,ul,  omega, 

1  w(19) , s (19) 

rear  fl,nl 

data  sf /I. 3 191, 1.9473,2 .5082,2.9529, 3 .2653, 3 .3462, 

1  3.0542,2.3640,1.5397, .6624, .21465, .73203,1.033, 

1  1.0222,  .68012, .27829,  . 3  7778 ,. 55739 ,  .46259/ 

data  ym/. 07246, .03741, .08306, .22726, .42925, .6717, 

1  .89061,  .99661,1.0366,  .9459,  .73952,  .47369,  .19894, 

1  .17588, .31972, .34974, .25216, .07516, .15072/ 

data  wee/. 293,  .361,  .433,  .509,  .588,  .679,  .756,  .845, 

1  .938,1.304,1.133,1.236,1.342,1.452,1.565,1.682, 

1  1.801,1.925,2.052/ 

data  wf/2 74, 277, 2 77, 2 80, 287, 29 0,293, 294, 292, 2 81, 163, 

1  122,110,101,87,40,313,290,280/ 

data  wm/169, 138,32, 19,22,23,23,26,27,24, 19, 10,339, 

1  242, 211, ].97, 185, 146,24/ 

C 

nodel=  33.27 
Im  =  5 . 

1  =  527.8 

rho  =  1.9905 

g  =  32.2 

ul  =26.0 
alpha=  8.1*. 001 
beta  =-  .74 
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ised  =  17 
omega=  2 . 0 
time  =  time*l/ul 
c 

c  Pierson-Moskowitz  Spectrxam 

c 

do  1  1  =1,19 

raof ( i) =sf ( i) *model*g/lm 
raom(i)  =yra(i)  *model*g 

wl (i) = (l+sqrt (l-4*ul*cos (2.6l8-psi)/g) )/ 

1  (2*ul*cos (2 . 618 -p&i) /g) 

w2 (i) = (1-sqrc (l-4*ul*cos (2.6l8-psi)/g))/ 

1  (2*ul*cos (2 . 618 -psi) /g) 

if  (wl (i) *w2 (i) .It .0)  then 
if  (wl(i).gt.O)  w(i)=wl(i) 
if  (w2(i).gt.O)  w(i)=w2(i) 
endif 

if  (wl (i) *w2 (i) .gt . 0)  then 
if  (abs (wl (i) ) .gt .abs (w2 (i) ) )  w(i) =abs (wl (i) ) 
if  (abs (w2 (i) ) .gt .abs (wl (i) ) )  w (i) =abs (w2 (i) ) 
endif 

s(i)  = (alpha* (g**2) / (w(i) **5) ) * 

1  exp (beta* (g/ (ws*w(i) )) **4) 

sx(i)  =s (i) / (1- (2*w(i) /g) *ul*cos (2 .618 -psi) ) 
vecx ( i) =sqrt (2*3x ( i) *raof (i) ) * 

1  cos (wee (i) *time-wf (i)*3. 14/180) 

vecm(i) =sqrt (2*sx(i) *raom(i) ) * 

1  cos (wee (i) *time-wm(i) *3 . 14/180) 

1  continue 

c 

c  Force  and  Moment  calculation 

c 

call  trap (19 , vecx, wee, fl) 
call  trap (19 ,vecm,wee,nl) 
fl=fl/ ( .5*rho*l**2*ul**2) 
nl=nl/ ( . 5*rho*l**3*ul**2 ) 
return 
end 
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